      program id_EU
      implicit real*8 (a-h,o-z)
      parameter(ndec=12,npat=4096,ngrid=51)
c
      dimension pa(ndec,2),pe(ndec,2),tl(npat+1,3)
      dimension bgrid(ngrid),g(ngrid),pz(ngrid)
      dimension er1(ngrid,ngrid),er2(ngrid,ngrid)
      dimension in(npat,ndec)
      common /pat/in
      write(*,'('' ***** id_EU *****'',20x,
     &/"    Ambiguity-Averse vs. EUM*")')
c
      open(unit=10,file='pat4096.data')
      do 10 n=1,npat
         read(10,*)idum,(in(n,k),k=1,ndec)
c         write(*,'(i5,x,12i2)')n,(in(n,k),k=1,ndec)
   10 continue
c
      gam = 2.d0*dlog(999.d0)
      call fpx(gam,1.d0,pe)
c
      do 20 j=1,ngrid
         if (j.eq.1) then
            pz(j) = 0.999d0
         else
            pz(j) = 0.5d0 + 0.01d0*(ngrid-j)
         endif
         g(j) = 2.d0*dlog(pz(j)/(1.d0-pz(j)))
         bgrid(j) = 0.d0 + 0.02d0*(j-1)
   20 continue
c
      do 100 j=1,ngrid
         gam = g(j)
         do 50 k=1,ngrid
            beta = bgrid(k)
            call fpx(gam,beta,pa)
c
            call fl(pa,pe,tl)
            call sort(npat+1,tl)
c
            cume = 1.d0
            cuma = 0.d0
            do 30 n0 = 1,npat+1
               if (tl(n0,1) .ge. 0.d0) goto 31
               cume = cume - tl(n0,3)  !cume are gt tl(n0,1)
               cuma = cuma + tl(n0,2)  !cuma are le tl(n0,1)
c               write(*,'(i4,5g13.4)')n0,(tl(n0,m),m=1,3),cuma,cume
   30       continue
   31       er1(j,k) = cuma             !cuma are le tl(n0-1,1) lt 0
            er2(j,k) = cume             !cume are ge 0
c               write(*,'(i4,5g13.4)')n0,(tl(n0,m),m=1,3),cuma,cume
            if(er1(j,k).gt.er2(j,k)) er2(j,k) = er1(j,k)
   50    continue
  100 continue
      er2(1,ngrid) = 0.5d0      ! when AA = EUM
c
      write(*,'(//"Max of Type I and Type II Errors:"/)')
      write(*,'(12x,51g12.4)')(bgrid(k),k=1,ngrid)
      do 110 j=1,ngrid
         write(*,'(52g12.4)')pz(j),(er2(j,k),k=1,ngrid)
  110 continue
c
      end
c***********************************************************************
      subroutine fpx(gam,beta,px)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12)
      dimension px(ndec,2)
c
      y1 = 10.d0
      y2 = 12.d0
      y3 = 15.d0
      y4 = 5.d0
      y5 = 6.d0
      y6 = 8.d0
      pjb = (1.d0+beta)/3.d0
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(1,1) = za/(za+zb)
      px(1,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(2,1) = za/(za+zb)
      px(2,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(3,1) = za/(za+zb)
      px(3,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(4,1) = za/(za+zb)
      px(4,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(5,1) = za/(za+zb)
      px(5,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(6,1) = za/(za+zb)
      px(6,2) = zb/(za+zb)
c
      za = y1*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(7,1) = za/(za+zb)
      px(7,2) = zb/(za+zb)
c
      za = y2*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(8,1) = za/(za+zb)
      px(8,2) = zb/(za+zb)
c
      za = y3*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(9,1) = za/(za+zb)
      px(9,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y4*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(10,1) = za/(za+zb)
      px(10,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y5*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(11,1) = za/(za+zb)
      px(11,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y6*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(12,1) = za/(za+zb)
      px(12,2) = zb/(za+zb)

c      do 10 k=1,ndec
c         write(*,'(i2,2g12.4)')k,(px(k,j),j=1,2)
c   10 continue
c
      return
      end
c***********************************************************************
      subroutine fl(pa,pe,tl)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12,npat=4096)
      dimension pa(ndec,2),pe(ndec,2),tl(npat+1,3)
      dimension in(npat,ndec)
      common /pat/in
c
      y1=0.d0
      y2=0.d0
      pt = 0.d0
      do 50 n=1,npat
         ta = 0.d0
         te = 0.d0
         do 10 k=1,ndec
            if (in(n,k) .eq. 1) then
               if (pa(k,1).gt.1.d-302) ta = ta + dlog(pa(k,1))
               if (pe(k,1).gt.1.d-302) te = te + dlog(pe(k,1))
            else
               if (pa(k,2).gt.1.d-302) ta = ta + dlog(pa(k,2))
               if (pe(k,2).gt.1.d-302) te = te + dlog(pe(k,2))
            endif
   10    continue
         tl(n,1) = ta-te
         tl(n,2) = dexp(ta)
         tl(n,3) = dexp(te)
         pt = pt + tl(n,3)
         y1 = y1 + ta*tl(n,2)
         y2 = y2 + te*tl(n,3)
c      write(*,'(i5,2x,5g12.4)')n,ta,te,tl(n,1),tl(n,2),tl(n,3)
   50 continue
      tl(npat+1,1) = 1.d2
      tl(npat+1,2) = 0.d0
      tl(npat+1,3) = 0.d0
c      write(*,'("pt = ",3g12.5)')pt,y1,y2
c      stop
c
      return
      end
c***********************************************************************
      subroutine sort(n,tl)
      implicit real*8 (a-h,o-z)
      dimension tl(n,3),a(n),b(n),c(n),indx(n)
c
      do 5 i=1,n
         a(i) = tl(i,1)
         b(i) = tl(i,2)
         c(i) = tl(i,3)
    5 continue
      call indexx(n,a,indx)
c
      do 10 i=1,n
         j = indx(i)
         tl(i,1) = a(j)
         tl(i,2) = b(j)
         tl(i,3) = c(j)
c      write(*,'(i5,x,3g12.4)')i,tl(i,1),tl(i,2),tl(i,3)
   10 continue
      return
      end

c***********************************************************************
      SUBROUTINE INDEXX(N,ARRIN,INDX)
      implicit real*8 (a-h,o-z)
      DIMENSION ARRIN(N),INDX(N)
      DO 11 J=1,N
        INDX(J)=J
11    CONTINUE
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          INDXT=INDX(L)
          Q=ARRIN(INDXT)
        ELSE
          INDXT=INDX(IR)
          Q=ARRIN(INDXT)
          INDX(IR)=INDX(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            INDX(1)=INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
          ENDIF
          IF(Q.LT.ARRIN(INDX(J)))THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        INDX(I)=INDXT
      GO TO 10
      END


